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Abstract:  Detecting  the  number  of  signals  and  estimating  the  parameters  of  the  sig¬ 
nals  is  an  important  problem  in  signal  processing.  Quite  a  number  of  papers  appeared 
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perturbation  analysis  of  the  data  auto  correlation  matrix  to  estimate  the  number  of 
sinusoids,  which  is  in  some  sense  a  subjective-based  method.  Recently  Reddy  and 
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H'lC  QUALITY  INSPECTED  S 


1  Introduction 


The  problem  of  detecting  the  number  of  terms  and  estimating  the  parameters  of 
sinuoids  in  presence  of  noise  is  an  important  problem  in  signal  processing.  In  the  last 
twenty  years  several  methods  have  been  proposed  to  estimate  the  frqucncics  of  the 
sinusoids  very  efficiently,  see  [1],  but  not  much  attention  has  been  paid  in  estimating 
the  number  of  sinusoids  presents  in  any  particular  signal.  Tufts  and  Kumaresan  [7] 
proposed  the  modified  forward  backward  linear  prediction  (MFBLP )  technique,  which 
is  capable  of  estimating  the  poorly  separated  frquencies,  with  short  data  lengts  and 
moderate  SNR,  but  unfortunately  it  is  known  to  give  inconsistent  estimates  (see  [10] 
and  [13]).  They  have  also  given  a  graphical  technique  to  estimate  the  number  of 
sinusoids,  which  is  very  subjective.  Some  of  the  other  techniques  which  are  available 
in  the  literature  for  example  [8],  [9],  [3]  can  be  used,  but  they  also  depend  on  the 
subjective  choice  of  the  individual  and  therefore  makes  it  difficult  to  implement  in 
practice. 

Rao  [10]  proposed  an  information  theoretic  criterion  to  estimate  the  number  of 
signals  for  undamped  exponential  model  but  it  is  observed  [5]  that  Rao’s  suggestion 
is  very  dufficult  to  implement  in  practice.  A  practical  implementation  procedure  has 
been  suggested  in  [5]  but  it  is  observed  that  the  proposed  method  depends  very  much 
on  the  penalty  function  used.  Some  suggestions  about  the  penalty  function  have 
been  given  mainly  based  on  the  computer  simulation  but  not  on  any  proper  analysis. 
Under  the  assumptions  of  normality  of  the  error  components,  recently  Reddy  and 
Biradar  [2]  proposed  two  criteria  one  is  AIC  type  and  the  other  one  is  MDL  type 
based  on  the  singular  value  decomposition  technique  and  using  both  the  forward  and 
backward  linear  predictor  equation.  They  also  obtained  the  probability  of  over  as 
well  as  under  estimation  following  the  assumption  made  by  [4]  and  [11].  They  use 
the  assumptions  of  normality  of  the  error  random  variables  in  developing  the  criteria 
as  well  as  in  the  performance  analysis.  The  performance  looks  quite  satisfactory. 

In  this  paper,  we  develop  a  method  mainly  based  on  the  extended  order  modelling 
and  singular  value  decomposition  technique.  We  use  penalty  functions  similarly  as 
AIC  and  MDL  but  instead  of  any  particular  penalty  function  a  class  of  penalty 
functions  satisfying  some  special  properties  has  been  used.  We  prove  that  any  penalty 
function  from  that  particular  class  will  give  consistent  estimate.  We  only  assume 
that  the  errors  are  independent  and  identically  distributed  (i.i.d.)  with  mean  zero 
and  finite  variance  to  prove  the  consistency  result.  In  fact  we  do  not  need  any 
distributioinal  assumptions  on  the  error  random  variables.  Next  we  obtain  an  estimate 
on  the  probability  of  wrong  detection  under  the  assumptions  that  the  errors  are 
i.i.d.  Gaussian  random  variables  with  mean  zero  and  finite  variance.  We  obtain  the 
probabilty  of  wrong  detection  using  matrix  perturbation  technique  and  large  sample 
approximation.  Once  we  obtain  the  probability  of  wrong  detection  for  any  particular 
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penalty  function,  we  choose  that  penalty  function  for  which  the  probability  of  wrong 
detection  is  minimum.  Some  simulations  are  performed  to  see  the  usefulness  of  our 
method  and  to  compare  it  with  the  other  existing  ones. 

The  rest  of  the  paper  is  organised  as  follows.  In  Section  2,  we  give  the  estimation 
procedure  and  the  consistency  result  is  proved  in  Section  3.  The  performance  analysis 
is  carried  out  in  Section  4  and  some  numerical  experiments  are  performed  in  Section 
5.  The  choice  of  the  penalty  function  is  suggested  in  Section  6  and  finally  we  draw 
the  conclusions  of  our  work  in  Section  7. 


2  Estimation  Procedure 

Consider  N  uniformly  spaced  data  samples  y{n)  of  M  real  sinusoids  corrupted  by 
additive  white  noise 

M 

y(n)  =  ^  aiSin(2'KUJin  +  (j>i )  +  ^(^)?  n  =  1, N  (1) 

i=l 


where  a^s  are  the  amplitudes,  can  be  any  arbitrary  real  numbers,  u's  are  the 
angular  frequencies  lying  between  0  and  .5,  0's  are  the  initial  phases  of  the  i  th 
frequency  lying  between  0  to  27T.  e(n)’s  are  i.i.d.  random  variables  with  mean  zero 
and  finite  variance  a2.  M,  the  number  of  signals,  is  assumed  to  be  unknown.  Given 
a  sample  of  size  N,  the  problem  is  to  estimate  M.  Let  us  assume  that  the  number 
of  sinusoids  can  be  at  most  K,  which  is  known,  i.e.  M  <  K.  Consider  the  following 
data  matrices; 

(  y{ 1)  •••  y(L)  \ 

Ap  =  :  :  :  ,  (2) 

v  y(N  -  L  +  l)  ...  y(N)  / 

and 


Rn  = 


(AT- L  +  l) 


[AyAp] 


here  L  is  any  integer  such  that  2 K  <  L  <  N  —  2 K  +  1. 
be  Ai  A L-  Compute 


(3) 


Let  the  eigenvalues  of  Rn 


IC(k)  =  log(A2*+i  +  1)  +  kCw-L+i 


(4) 
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for  k  =  1,  . . . ,  2 K,  here  C n  (penalty  functions)  satisfies  the  following  conditions 

(a)Cn  >  0  (6)C„^0  («)Jgg,-M»  (5) 

Then  our  criterion  is  the  following;  Choose  M  such  that 

IC{M)  =  rnin[IC(  1), ....  IC(K)}  (0) 

So  M  is  an  estimate  of  M.  In  the  next  section  we  prove  that  M  converges  to  M 
almost  surely  for  any  CV,  satisfying  (5),  under  the  assumptions  that  the  errors  are 
i.i.d.  with  mean  0  and  finite  varinace  o 2  >  0. 


3  Consistency  Results 


Let  us  denote  the  (p,q) th  element  of  the  matrix  AjAp  by  apq.  Therefore 


apq 


N-L 


53  y(p + s)2/(? + 

a= 0 


(7) 


Now  write; 

» 

M 

y(n)  =  ^2  aisinfawin  +  <f>i)  +  c(n) 

t=i 

M 

=  ^2  -r[exp(j27ra;in)exp(j^i)  -  exp(-j2iruin)exp(-j(t)i)}  +  e(n) 
i=x  *3 

2M 

=  ^2  Aiexp(j2Trqin)  +  e(n)  (8) 

i=i 

where  j  =  y/— 1,  Ai  =  ^  e^1,  A  2  =  —  -y*  e  =  ^1,  P2  —  “^i)  an<^  the  others 

are  defined  similarly.  Since  y{q  +  s)  =  y{q  4-  s )  (*-’  denotes  the  complex  conjugate), 

AT-L 

aP9  =  53  y(p+s)y(<i  +  s) 

3=0 

AT-L  2M 

=  EE  Aiexp(j2icr]i(p  +  s))  +  e(p  +  s))  x 

s=0  i=l 
2Af 

(J3  Aiexp{-j2'KT)i{q  +  s))  +  e(q  +  s)) 

i= 1 

N-L  2 M  2 M 

=  E(E  Aiexp(j2'Kr]i(p  +  s)))(53  Aiexp(—j2irr)i(q  +  s))) 

s=0  t=l  «=i 
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N-L 


2M 


+  X  e(p  +  s)(£  Aiexp{-j2-KT]i{q  +  5))) 

«=o  t=l 

2M 

+  2  c(9  +  s)(X  Aiexp(j2irqi{p  +  s))) 


1=1 


5=0 

N—L 

+  XI  eb  +  s)«(g  +  s) 

s= 0 


(9) 


vve 


Therefore  by  the  law  of  the  iterated  logarithm  of  the  M  -  dependent  sequence, 
can  say  that  for  fixed  L  as  N  tends  to  infinity,  we  have 

R"  =  TA?Ar  =  +  A  +  0(M°g(^  +  1))'/2(^.)  (10) 

here  ‘a.s.’  stands  for  almost  surely  and  the  matrix  ft  and  D  are  as  follows: 

(  exp(27T7?1)  ...  exp(2'KT]2M)  \ 

, D  =  diag(||  yl1||2)...,||  A2Mf)  (11) 


ft  = 


V  exp(2nr]lL)  . . .  exp(2nr)2M L) 


Let  us  denote  the  matrix  ftDft^  by  E  and  let  the  non  zero  eigenvalues  of  E  be 
Pi  >  ...>  faM-  Let  us  denote  the  ordered  eigenvalues  of  E  +  cr2I  by  Ax  >  . . . \2M  > 
X™+i  =  ...  =  XL  =  a2,  A i  =  pi  +  a2  for  i  =  1, . . . ,2M  and  A.  =  a2  for  i  = 
2M  +  1, . . . ,  L.  We  need  the  following  result  for  further  development: 

Lemma  1  :  Let  P  and  Q  be  two  Hermitian  matrices  of  order  m  x  m  and  let  the 
spectral  decomposition  of  P  and  Q  be  as  follows 


m  m 

P  =  X5iPiPi*  and  Q  =  X  TiQiqf1  (12) 

i=l  i=l 

where  <  ...  <  5m  and  71  <  ...  <  7m  are  the  ordered  eigenvalues  and  pj’s, 
qx’s  are  the  corresponding  orthonormal  eigenvectors  of  P  and  Q  respectively,  then  if 

bi*  —  Qik\  <  ct  for  all  i,  k  =  1, . . . ,  m  and  for  some  a,  then  there  is  a  constant  C  such 
that  ft  -  7i|  <  Cot. 

Proof  :  It  mainly  follows  from  von-Neumann’s  inequality  but  see  [12]  for  details. 
Using  Lemma  1  and  (10)  we  obtain 


X  _  X  .  „,loglog(N-L  +  1),1/2,  x 
A,-Ai  +  0(-  jv_L+1 - )  (a.,.) 


(13) 


Now  to  prove  M  is  a  consistent  estimate  of  M  observe  that  it  is  enough  to  prove 
IC(q,  Cn-.l+i)-IC(M,  CN-L+i)  >  0;  for  q  —  1, ... M-l, M+l, ... K  (14) 
as  N  — >  00. 
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Case  I:  q  <  M 

IC(q,  Cn-l+i)  —  IC(M,  Cm-l+i) 


Case  II:  q  >  M 

IC(q,  CN.L+1)  -  IC(M,  CV-l+1) 


here  h,  =  =  h2. 

Now  observe  that; 


l°g(^2g+l)  —  log(A2M+l)  +  (<7  ~  M)  C^-L+l 
l°g(^2<j+l  +  1)  —  log(A2M+l  +  1) 
log(/i2M+i  +  a2  +  1)  -  log(cr2  +  1)  >  0  (15) 


l°g(^29+l)  —  log(A2M+l)  +  (q  —  M )  Ctf-L+1 
log(l  +  cr2  +  hi)  -  log(l  +  a2  +  h2) 


+(q  —  M )  Cn-l+ i 


log(l  +  cr2)  +  hi 


1 


1  +  cr 


—  h2 


1  +  a2 


-  log(l  +  cr2)  +  {q  -  M)  Cn-l+i 
/o^off(iV  -  L  +  1) 

1  N-L  +  l  ’ 


(16) 


(IC(q,  Cn-l+i)  -  IC(M,  CN.L+l)) 


=  (q~M)  + 

1  ~  flog  log(iV  —  L  +  1) 


Cn-L+i 

1 


o(- 


N-L  +  l 


)l/2 


Cn- 


L-f  1 


oi1*^;11)  in) 


From  the  properties  of  Cn  (see  (5)),  it  follows  that  the  second  and  the  third  term  of 
the  right  hand  side  of  (17)  goes  to  zero  as  N  tends  to  infinity,  therefore  (17)  implies 

TT - - (IC(q,  Cn-l+i)  -  IC(M,  CN.L+i))  -+  (q  -  M)  >  0  (18) 

^N-L+ 1 


Therefore  we  can  conclude  that  for  large  (N  —  L  4- 1) 


(IC(q,  CN.L+i)  -  IC(M,  CN.L+i))  >  0  (19) 


Combining  (15)  and  (19)  we  obtain  (14). 


4  Performance  Analysis 


In  this  section  we  obtain  an  upper  bound  for  P(M  ^  M).  Now 
P{M^M)  =  P(M  <  M)  +  P(M  >  M) 


6 


(20) 


M-l  ~  K 

=  £  P{M  =  q)  +  P{M  =  q) 

9—0  q=M+\ 

A/-1 

=  E  p(ic(q,  cN_L+1)  -  ic(m ,  cN_L+1)  <  o) 

9=0 

K 

+  E  ^(^(9>  CN.L+1)-IC(M,  cN_L+l)<0) 

q=M+ 1  '  ’ 

Let  s  consider  two  different  cases: 

Case  I,  q  <  M 

P{IC(q,  CN-L+i)  -  IC(M,  CN_L+l)  <  0) 

=  P(log(X2q+i  +  1)  -  log( X2M+1  +  1)  +  (q  -  M)  CN-L+l  <  0) 

=  P(log(X2q+l  +  1)  -  log( X2M+1  +  l)  +  (q-M)  CN.L+1  < 

log{X2M+1  +  1)  -  log(X2M+1  +  1)  +  log(X2q+l  +  1)  -  log(X2q+1  +  1)) 

<  P(log{ X2q+1  +  1)  -  log{ X2M+l  +  1)  <  (M  -  q)  CN.L+l  + 

\log(X2M+i  +  1)  -  log{ X2M+l  +  1)|  +  \log(X2q+1  + 1)  -  log(X2q+l  +  1)|)  (21) 

Let  6  be  such  that  for  large  N 

l°9(^2q+i  +  1)  -  log(X2M+i  +  1)  >  (M  -  q)  CN-L+1  +  5  (22) 

Therefore  for  large  N 

P(IC(q,  Cn_l+ i)  —  IC(M,  Cn-l+x)  <  0)  =  0  (23) 

Case  II,  q  >  M 

P{IC{q,  Cn-l+i)  ~  IC(M ,  Cn-l+i)  <  0) 

—  P(log(A29+1  +  1)  —  log(A2M+i)  +  (q  —  M )  Cat_£,+i  <  0) 

=  -P(log(A2M+1  +  1)  -  log(A2?+1  +  1)  >  (q  -  M)  CN_L+1)  (24) 

Therefore  observe  that  from  (23)  it  is  clear  that  for  large  N  the  probability  of  under¬ 
estimation  is  zero  and  to  obtain  (24)  we  need  to  know  the  joint  distributions  A*  for 
k  =  2M  + 1, . . . ,  L.  Let’s  assume  at  this  point  that  e(n)’s  are  i.i.d.  Gaussian  random 
variables  with  mean  zero  and  variance  a2.  Later  on  we  see  that  it  is  possible  to  relax 
this  assumption.  Let’s  denote  the  matrix  (E  +  a2I)  by  R.  Therefore  from  the  central 
limit  theorem,  we  can  say  that  asymptotically 

ViV  -L  +  l(Hec(RN_L+1)  -  Kec(R))  (25) 

will  be  normally  distributed  with  mean  vector  zero  and  certain  L2  x  L2  variance 
covariance  matrix  T  (say).  Here  Vec(.)  of  an  Lx  L  matrix  denotes  the  L2  x  1  vector 
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stacking  the  columns  one  below  the  other.  Using  the  perturbation  theory  ([15];  page 
66)  let’s  write 

Rn-l+i  =  R  +  (Rn-l+i  —  R)  =  R  +  eB  (26) 

Here  0  <  £  <  1,  and  B  is  a  Hermitian,  zero  mean  matrix  with  elements  that  are 
asymptotically  jointly  Gaussian,  which  follows  from  (25).  Let  \  be  any  particular 
eigenvalue  of  R  and  let  the  corresponding  normalized  eigenvector  be  Zj.  Observe 
that  if  A i  is  a  multiple  eigenvalue  of  R,  then  z\  is  not  unique,  but  then  we  take  one 
particular  Zi  (see  [15]  page  69).  Let  A  *  be  the  corresponding  perturbed  eigenvalue  of 
Rn-l+i,  then  from  ([15]  page  69)  it  follows  that 

Aj  =  Ai  +  e(zjrBzj)  (27) 

Since  the  elements  of  B  are  normally  distributed,  therefore  z^Bzi,  which  is  a  linear 
combination  of  the  elements  of  B,  will  also  be  normally  distributed  with  mean  zero 
and  finite  varinace.  Now  from  (27)  it  follows  that 

E(\i)  ~  Aj  (28) 

Now  to  compute  the  covariance  between  Aj  and  Aj,  let  z;  and  zj  be  any  two  orthonor¬ 
mal  eigenvectors  corresponding  to  Aj  and  Aj  respectively,  therefore  we  have 

Aj  =  Aj  +  e(zjBzi)  (29) 

Aj  =  Aj  +  e(zfBzj)  (30) 

Therefore 

E(\i  -  Aj)(Aj  -  Aj)  =  c2£'(z?’Bzi)(zjrBzj)  (31) 

Since  we  are  mainly  interested  about  the  distribution  of  the  repeated  eigenvalues,  so 
let’s  take  z\  and  Zj  be  any  two  orthonormal  eigenvectors  corresponding  to  a2,  therefore 
we  can  obtain  after  some  simplifications; 


E  (zfBziKzjBz,) 


(N  —  L  +  1)V 


N-L+l  JV—L*f  1 

E  £  (zf r PqZj )  (zJrqpZi) 


P=  1  9=1 


(32) 


Here  rp9’s  are  Lx  L  matrices  and  they  are  defiend  as  follows  TP9  =  T^p  and  rp9  =  0 
if  ( q-p )  >  L.  For  q  >  p, Tpq  has  all  the  elements  to  be  zero  except  at  the  positions 
(q  —  p  +  1,1), . . .  ,(L,p  —  q  +  L)  which  are  ones,  as  follows; 


/  0  ...  0  0  \ 


r 


P9  ~ 


1  ...  0  0 
01  ...  0  ’ 


\o  ...  1  0/ 


(33) 
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Therefore  we  obtain  that  A2M+i  and  A2,+i  for  q  >  M  are  jointly  normal  each  of  them 
with  mean  a2  and  it  has  the  variance  covarinace  matrix  which  can  be  obtained  horn 

(32). 

Observe  that  even  if  we  know  the  joint  distribution  of  K  for  i  =  2M  +  1,  •  •  • ,  L, 
theoretically  it  is  very  difficult  to  compute  (24).  We  use  simulation  technique  similarly 
as  [11]  to  compute  (24)  by  using  the  joint  distribution  A2M+i  and  A2(?+1.  the  details 
will  be  explained  in  the  next  section. 

Observe  that  since  we  assume  that  the  eigenvalues  of  E  are  distinct  therefore 
the  distribution  of  Ax, ... ,  A2M  will  be  independent  of  A2M+i,  •  •  • .  xl  asymptotica  y. 
Moreover  they  will  be  jointly  normal  with  mean  Ai  and  variance  ^lTT  for  z  “ 
1, . . . ,  2 M  see  [14].  and  they  will  be  independent  of  each  other. 


5  Numerical  Experiments 


In  this  section  we  perform  some  numerical  experiments  to  present  both  the  effective- 
necs  of  cm  method  and  the  usefulness  of  the  analysis.  We  consider  the  same  mode. 

as  that  of  [2]. 


The  data  sample  are  generated  from  the  following  model: 
y{n)  =  y/20sin(2nwin)  +  \/20sm(27ru;2n  +  <f>)  +  e(n); 


n  =  1, . . . ,  N. 


here  o>,  =  .2  and  <v2  =  .2  +  6  with  <5  =  j,  and  N  =  64.  Here  e(n)'s  are  i.i.d.  Gaussian 
random  variables  with  mean  zero  and  variance  o2,  which  “  yWe 

eive  the  desired  signal  to  noise  ratio  (SNR)  defined  as  SNR  lOlogw  /  • 

use  twelve  different  C„,  all  of  them  satisfying  (5)  but  converging  to  zero  a  different 

rates.  We  define  them  as  C$> •  They ;  are  as  follows;  CN  =  (77)  ,C„- 
(^)‘2>Ci?  =  (f)'3’Cn  =  (^)'4-Cn  =  to{j ~  (iogiv)'2>C^i2)  N 

take  K  =  8,  i.e.  the  maximum  number  of  sinusoids,  same :  as  ]  e 

Out  of  500  simulations,  the  percentage  of  correct  estimate  (PCE),  th P*  c  g^ 
of  under  estimate(PUE)  and  the  percentage  of  over  estimate  (POE)  a  P 
for  SNR  =  5dB  and  SNR  =  lOdB.  We  also  obtain  the  theoretical  value  for 
upper  bound  of  the  probability  of  over  estimate  as  follows.  We .draw r  a  samp 
size  (L-2M)  from  Gaussian  random  variable  with  mean  a  and  variance  covananc 
mate  given  by  (32).  we  order  them  as  A.  for  i  =  2 M  L  and  Id eck  whetlm 
n  \  -  Innln  A  >  (a  -  M)  Cn-l+1  is  true  or  not  for  q  =  M  +  1,  •  •  • . 

Lpeatthis  process^OOO  times  and  compute  the  percentage  of  times  it  is  true  and t 
gives  an  estimate  of  (24).  Observe  that  although  we  have  seen  that  the  pro  a  1 1  y 
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underestimate  will  be  zero  for  large  sample  size,  we  have  calculated  the  probability 
of  underestimate  using  the  distribution  properties  of  the  \u  . . . ,  A2M+1  and  using  the 
same  procedure  as  above  by  repeating  over  5000  times.  The  results  are  reported  in 
Table  1  for  SNR  =  5dB  and  in  Table  2  for  SNR  =  lOdB.  The  quantity  within  the 
bracket  indicate  the  theoretical  bounds  of  POE’s  and  PUE’s. 


Table  1  Table  2 


PUE 

PCE 

POE 

cil' 

0(0) 

97 

3(5) 

421 

0(0) 

84 

16(22) 

4S) 

0(0) 

59 

41(50) 

44) 

0(0) 

31 

69(82) 

UAT 

0(0) 

49 

51(56) 

0(0) 

100 

0(3) 

<#> 

0(0) 

94 

6(14) 

ct8) 

0(0) 

82 

18(23) 

cf 

0(0) 

64 

36(39) 

cr 

0(0) 

94 

6(12) 

eft1’ 

0(0) 

29 

71(82) 

r(12) 

0(0) 

98 

2(6) 

At  finite  sample  size  the  performance  of  the  proposed  method  very  much  depends 
on  the  penalty  function  used,  although  all  of  them  give  consistent  estimates  as  the 
sample  size  tends  to  infinity.  From  Table  1  and  Table  2  it  is  clear  that  the  performance 
of  all  the  methods  becomes  worse  at  low  SNR,  which  is  not  very  surprising.  It  is 
important  to  observe  that  the  theoretical  probability  matches  quite  well  in  almost 
all  cases  considered  and  the  estimates  are  better  in  most  of  the  cases  at  high  SNR. 
It  seems  that  the  approximations  made  for  the  large  sample  work  reasonably  well  at 
moderate  sample  sizes  and  for  mederate  SNR. 


6  How  to  Choose  the  Penalty  Function? 


Looking  at  the  tables  it  is  clear  that  the  theoretical  bounds  are  quite  close  to  the 
actual  one.  But  unfortunately  without  knowing  the  original  parameters  we  can  not 
calculate  the  theoretical  probabilites.  Reddy  and  Biradar  [2]  also  did  not  raise  this 
question  that  how  to  estimate  these  bounds.  One  way  of  course  by  replacing  these 
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parameters  values  by  their  estimates.  We  would  like  to  estimate  these  probabilities 
with  the  help  of  the  given  sample  and  use  it  to  choose  the  proper  penalty  function. 

e  idea  is  as  follows:  From  any  particular  realization  of  the  model,  we  compute 
ne  matrix  Rn-l+i  and  obtain  the  eigenvalues  and  the  corresponding  eigenvectors. 
Now  suppose  using  the  penalty  function  C$\  we  estimate  the  order  of  the  model 
as  M„.  Assuming  M„  is  the  correct  order  model  we  compute  the  estimate  of  o2 
byaveraging  the  last  L  -  2  Mt  eigenvalues.  Now  using  the  joint  distribution  A,  for 

M.k.+  7 1 '  •  ’ ■L’  (24)  can  be  calculated  by  using  the  simulation  technique  as 
discussed  in  the  previous  section.  Similarly  from  A,, . . . ,  A2„,+1,  we  can  compute 

an  estimate  of  the  probability  of  under  estimate.  Therefore  adding  the  two,  we  can 
o  tain  an  estimate  of  the  upper  bound  of  the  probability  of  wrong  detection.  It  can 
be  shown  easily  that  for  large  N,  the  estimate  of  the  probabiltiy  of  wrong  detection 
under  the  asuumption  of  correct  order  model  will  be  less  than  the  estimate  of  wrong 
detection  under  the  assumption  of  lower/higher  order  model,  because  the  former  one 
goes  to  zero  as  N  tends  to  infinity,  where  as  the  later  one  goes  to  a  positive  quantity. 


We  use  this  idea  and  compute  the  estimate  of  the  probabiltiy  of  wrong  detection 
for  all  the  criteria  and  choose  that  one  which  gives  the  lowest  estimate  of  the  proba¬ 
bility  of  wrong  detection.  We  have  used  the  same  model  and  the  same  set  of  penalty 
functions  and  in  each  trial  we  choose  that  penalty  function  which  gives  the  lowest 
estimate  of  the  probabilty  of  wrong  detection.  In  each  trial  we  run  100  simulation 
to  compute  the  estimate  of  the  probability  of  error  and  it  is  repeated  over  500  trials. 
The  result  is  reported  in  Table  3,  which  indicates  the  percentage  of  underestimate, 
correct  estimate  and  over  estimate  out  of  these  500  replications. 


Table  3 


SNR 

PUE 

PCE 

POE 

5dB 

0 

94 

6 

lOdB 

0 

95 

5 

From  Table  3,  it  is  observed  that  the  proposed  method  works  quite  well.  For  the 
same  model  it  is  observed  [2]  that  when  SNR  =  5dB,  the  simulation  shows  that  MDL 
criterion  can  detect  at  most  88  percent  and  at  SNR  =  10dB  it  can  detect  at  most 
91  percent  correctly  but  in  our  case  it  is  observed  that  it  can  detect  94  percent  and 
95  percent  respectively.  It  is  higher  than  what  they  have  obtained,  which  is  not  very 
surprising  because  we  have  a  class  of  penalty  functions  and  from  there  we  are  trying 
to  choose  the  best. 

Another  question  naturally  comes  that  how  to  choose  the  class  of  penalty  func¬ 
tion.  We  suggest  the  following  way;  take  any  particular  class  of  resonable  size  may  be 
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around  ten  or  twelve  then  obtain  the  estimate  of  the  probability  of  wrong  detection 
for  all  the  methods  and  compute  the  minimum  value,  if  the  minimum  itself  is  high 
that  means  the  error  has  been  calculated  using  wrong  order  model  and  the  class  of 
penalty  function  is  not  good  and  on  the  other  hand  if  it  is  low,  that  means  this  class 
is  fine.  We  have  computed  the  mean  and  standard  deviation  of  the  minimum  proba- 
biltiy  of  error  when  the  model  has  been  chosen  correctly  and  when  the  model  has  not 
been  chosen  correctly  for  the  same  experiments.  The  results  are  reported  below  in 
Table  4a  and  Table  4b,  where  Table  4a  gives  the  mean  and  standard  deviation  (s.d) 
of  the  minimum  probability  of  error  when  the  model  has  been  chosen  correctly  and 
Table  4b  gives  the  result  when  the  model  has  not  been  chosen  correctly. 


Table4a 


Table4b 


SNR 

mean 

s.d 

5dB 

lOdB 

.2850 

.2540 

.1360 

.1288 

SNR 

mean 

s.d 

5dB 

lOdB 

.0489 

.0255 

.0845 

.0527 

From  the  Tables  4a  and  4b  it  is  clear  that  the  mean  of  the  minimum  probability  of 
error  under  the  assumptions  of  the  correct  order  model  is  much  smaller  than  under 
the  assumption  of  wrong  order  model.  So  it  is  expected  that  the  minimum  probabiltiy 
of  error  can  give  us  some  idea  whether  the  class  is  sufficient  or  not. 


7  Conclusions 

In  this  paper,  we  propose  a  method  to  detect  the  number  of  sinusoids  present  in 
a  sinusoidal  signal  corrupted  by  additive  noise.  We  use  extended  order  modelling 
and  singular  value  decomposition  technique.  We  prove  the  strong  consistency  of  the 
proposed  method  under  the  assumptions  that  the  errors  are  i.i.d.  with  mean  zero  and 
finite  variance.  We  have  done  some  performance  analysis  of  the  proposed  method 
under  the  assumption  that  the  errors  are  Gaussian.  But  theoretically  speaking  it  is 
possible  to  relax  this  condition.  Observe  that  (25)  is  true  asymptotically  even  if  the 
errors  are  not  Gaussian.  But  in  that  case  more  sample  size  is  required.  Another 
important  point  is  to  observe  that  we  can  use  both  the  backward  and  forward  data 
to  compute  RN  of  Section  2  and  in  that  case  it  is  not  very  difficult  to  prove  that  the 
proposed  criteria  will  be  consistent  but  unfortunately  in  that  case  the  performance 
analysis  becomes  difficult,  particularly  the  expression  (32)  will  have  more  terms  so  it 
is  not  pursued  here.  It  is  observed  that  our  method  behaves  quite  satisfactory  and 
works  better  than  the  existing  methods. 
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